############################################
############################################
# previous: 6-VAP_FutProj ##################
############################################

###### Summary ############################################################################################
#This script loads the future projections for habitat suitability for each modelling unit and summarises them the median, 
#90&, and 10% quantile across all scenarios, i.e. it calculates the median predicted future habitat suitability for each modelling unit
#as well as the highest (90%) and lowest (10%) expected future habitat suitability. For each, current, future median, future Q90 and 
#future Q10 rasters, it crops them to an area 10decimal degrees around known occurrences & aves them, then thresholds them & saves them.
#For future median, Q10, and Q90, it also calculates the change in suitability compared to current conditions & saves it. Lastly, the sfipt 
#creates .png images of maps of current & future habitat suitability and change therein.
###########################################################################################################

############################################
# next: 8-VAP_CDCut ########################
############################################

#######################################################################################################
##### MaxEnt3 Plotting ################################################################################
#######################################################################################################

##########
#Settings#
##########

#load packages:
library(R.utils)
library(raster)
library(maptools)
library(rgdal)
library(stringr)
library(sf)

#basic model settings
#set project name:
projectName<-"WHOSnakes"

#define variable selection metrics:
crits<-c("permutation.importance")   #c("permutation.importance","contribution","Test.gain.without","Test.gain.with.only","AUC.without","AUC.with.only")
critThresh<-c(1)                     #c(1,1,1,10,1,75)

#define directories
MyDir<-paste("/home/jc217070/Maxent",sep="")
SWDDir<-paste(MyDir, "/Maxent_SWDs/spp_swds",sep="")
OCCDir<-paste(MyDir, "/Maxent_SWDs/spp_occs",sep="")
OutDirCrossval<-paste(MyDir,"/Outputs2_Crossval",sep="")
OutDirFinal<-paste(MyDir,"/Outputs3_Final",sep="")
PlotDir<-paste(MyDir,"/Maxent_Data/Plots/Plots_WHOSnakes_2020_09_01",sep="")
dir.create(path=file.path(MyDir,"tmp.pbs"),showWarnings = FALSE)
PbsDir <- paste(MyDir, "/tmp.pbs", sep="")
setwd(PbsDir)

#read in lookup tables
DatSum<-read.table(file=file.path(paste(MyDir, "/Maxent_LUTables/",projectName,".csv",sep="")), header=TRUE, sep=",") #load data summary for all species
spp<-DatSum$gen_sp_subtax

mapspp<-unique(DatSum$ModUnit)
mapspp<-paste("ModUnit_", gsub(" ","_",mapspp), sep="")

#mapsppB<-unique(DatSum$ModUnit[DatSum$Redo==1])
#mapsppB<- (c("Dendroaspis polylepis","Echis ocellatus", "Echis leucogaster","Naja senegalensis","Naja savannula","Bitis arietans"))
mapsppB<-unique(DatSum$ModUnit)
mapsppB<-paste("ModUnit_", gsub(" ","_",mapsppB), sep="")

FutLayersDir<-paste("/home/jc217070/GIS_and_raw_Layers/final_layers/WorldClim2.1/future")
scenarios<-list.dirs(path=paste(FutLayersDir),full.names=FALSE,recursive=FALSE)
scenarios<- unique(gsub("_2041-2060","",gsub("_2081-2100","",scenarios)))
years<-c("2050","2090")
#years<-c("2050")

#read in base layers
worldmap <- st_read(paste(PlotDir,"/worldmap/WHO_countries.shp",sep=""))
WHOpolies <- st_read(paste(MyDir,"/WHOSnakesShapes/SSN_CTY_CAT_BRG_FEA_IN_Mar2023.shp",sep=""))
mybg<-raster(paste(PlotDir,"/Cop_FAPAR_mean_2010-2019_WW1km.asc",sep=""), crs="+init=epsg:4326")

#define plot colors
my_col1 = rev(terrain.colors(n = 20, alpha=1))
my_col2 = rev(terrain.colors(n = 20, alpha=0.5)) #alpha=0:1, higher = less transparent
my_col3 = c("white","green")
my_col4 = c("white","grey", "green")
my_col5<-c("royalblue4","royalblue3","royalblue2","royalblue1","steelblue1","lightcyan",
           "wheat1",
           "orange","darkorange","orangered","red","red3","red4")
mycols1<-c("grey90","grey85","grey80","grey75","grey70",
           "grey65","grey60","grey55","grey50","grey45",
           "grey40","grey35","grey30")
mycols2<-c("royalblue","cyan","navy","skyblue","gold",
           "red","darkred","darkorange","orangered",
           "deeppink","pink","purple","slateblue1","lightsalmon","burlywood4","coral3",
           "deeppink4","darkorange4","cyan4","purple4")
WHOcol <- rgb(150, 180, 250, max = 255, alpha = 0) #alpha=0 means completely transparent

##############
# main script#
##############


for (sp in mapspp) { #for each modelling unit...
  
  if (sp %in% mapsppB){

    n<-which(mapspp==sp) #gets row number for current species (for reference back to full data frame)
    spX<-gsub("ModUnit_","",sp)
    is<-which(DatSum$ModUnit==gsub("_"," ",spX))
    
    for (crit in crits){
      
      for (year in years){
        yearLong<-ifelse(year=="2050","2041-2060","2081-2100")

      
      # if (file.exists(paste(OutDirFinal, "/", crit,"/", sp, "/", sp, "_UKESM1-0-LL_ssp585_",yearLong,".asc.gz", sep="")) |
      #     file.exists(paste(OutDirFinal, "/", crit,"/", sp, "/", sp, "_UKESM1-0-LL_ssp585_",yearLong,".asc", sep=""))){ #safety, skip species with no output
      #   
        SpOutDirFin<-paste(OutDirFinal, "/", crit,"/",sp,"/future", sep="")
        
        if(!file.exists(paste(SpOutDirFin,"/",sp,"_",year,"_MedianCropThresh.asc.gz",sep="")) |
           !file.exists(paste(SpOutDirFin,"/",sp,"_",year,"_Q90CropThresh.asc.gz",sep=""))    |
           !file.exists(paste(SpOutDirFin,"/",sp,"_",year,"_Q10CropThresh.asc.gz",sep=""))    ){
          
        #get species records used for model:
        sprecs<-read.table(paste(OCCDir,"/",spX,"_mod.csv",sep=""), sep=",", header=TRUE)   #read in species recs
        recnum<-length(sprecs[,1])
        
        # read in the maxent results thresholds
        results <- read.csv(paste(OutDirFinal, "/", crit, "/", sp, "/maxentResults.csv", sep=""))
        
        # thresholds
        Oms<-read.csv(paste(OutDirFinal, "/", crit, "/", sp, "/",gsub("ModUnit_","",sp),"_omission.csv", sep=""))
        
        Omsn<-which(Oms$Training.omission==max(Oms$Training.omission[Oms$Training.omission<=0.05]))
        myt<-Oms$Corresponding.Cloglog.value[Omsn]
        
        # thresholds
        thresh1<-as.numeric(results[1,c("Maximum.training.sensitivity.plus.specificity.Cloglog.threshold")])
        thresh2<-as.numeric(results[1,c("Equate.entropy.of.thresholded.and.original.distributions.Cloglog.threshold")])
        thresh3<-as.numeric(results[1,c("Balance.training.omission..predicted.area.and.threshold.value.Cloglog.threshold")])
        thresh4<-myt
        
        threshs<-c(thresh1,thresh2,thresh3,thresh4)
        t <- 3
        #t<-as.numeric(DatSum$thresh[is])[1] #only use first one as they are all the same for each species within this mapping unit
        thresh<-threshs[t]

        
        if(length(sprecs[,1])>0){
          
          Sys.sleep(0.1)
          print(paste(Sys.time(), "making extent for modelling unit",sp, sep=" "))
          
          
          rangewidthDG<-max((max(na.omit(sprecs$long))-min(na.omit(sprecs$long))),
                            (max(na.omit(sprecs$lat))-min(na.omit(sprecs$lat))))
          rangewidthM<-rangewidthDG*100000 #turn dec deg into meters
          rangewidthM<-ifelse(rangewidthM>3000000,3000000,rangewidthM) #make buffer no bigger than this
          rangewidthM<-ifelse(rangewidthM<1000000,1000000,rangewidthM) #make buffer no smaller than this
          y<-rangewidthM
          
          #make extent object describing a box around occurrence records (margin= 4 times x and y extent of occurrence records):
          xrange<- (max(na.omit(sprecs$long)))-(min(na.omit(sprecs$long)))
          yrange<- (max(na.omit(sprecs$lat)))-(min(na.omit(sprecs$lat)))
          xcenter<- (min(na.omit(sprecs$long)))+(xrange/2)
          ycenter<- (min(na.omit(sprecs$lat)))+(yrange/2)
          rangemax<-max(xrange,yrange)
          xmin<-ifelse((xcenter-(rangemax/2)- y/100000) >= -180, (xcenter-(rangemax/2)- y/100000), -180)
          xmax<-ifelse((xcenter+(rangemax/2)+ y/100000) <=  180, (xcenter+(rangemax/2)+ y/100000),  180)
          ymin<-ifelse((ycenter-(rangemax/2)- y/100000) >=  -90, (ycenter-(rangemax/2)- y/100000),  -90)
          ymax<-ifelse((ycenter+(rangemax/2)+ y/100000) <=   90, (ycenter+(rangemax/2)+ y/100000),   90)
          myext<-extent(xmin,xmax,ymin, ymax)
        } else{
          myext<-extent(-180,180,-90,90)
        }
        
        
        #######################################
        ############### Summarizing ###########
        #######################################
        
        print(paste(Sys.time(), "calculating Med, Q90 & Q10 for", sp, sep=" "))
        
        #get file names for all future scenario ascs
        outfiles<-paste(SpOutDirFin,"/",sp, "_", scenarios,"_",yearLong,".asc",sep="")
        
        ascstack<-stack()
        for (file in outfiles){
          
          
          
          
          if(file.exists(file) |
             file.exists(paste(file,".gz",sep=""))){
            
          
          
          
          i<-which(outfiles==file)
          Sys.sleep(0.1)
          print(paste(Sys.time(), "stacking projection for scenario", scenarios[i], sep=" "))
          
          if(file.exists(paste(file,".gz",sep="")) & !(file.exists(paste(file)))){
          gunzip(filename=paste(file,".gz",sep=""), #unzip model if only zipped file exists
                 destname=paste(file,sep=""),
                 remove=FALSE, overwrite=TRUE)
          }
          
          layer<-raster(file)
          layer<-crop(layer, myext)
          ascstack <- stack(ascstack, layer)
          
          # if (file.exists(file)){
          # #zip future projection file after I'm done with it
          # Sys.sleep(0.1)
          # print(paste(Sys.time(),"zipping up raw future asci",sep=" "))
          # gzip(filename=paste(file),overwrite=TRUE, remove=TRUE)
          # }
          #   
        }}
        
        
        #threshold current SDM
        
        
        Sys.sleep(0.1)
        print(paste(Sys.time(), "Current", sep=" "))
        if(!(file.exists(paste(OutDirFinal, "/", crit,"/", sp, "/", sp, "_current.asc", sep="")))){
          gunzip(filename=paste(OutDirFinal, "/", crit,"/", sp, "/", sp, "_current.asc.gz", sep=""), 
                 destname=paste(OutDirFinal, "/", crit,"/", sp, "/", sp, "_current.asc", sep=""),
                 remove=FALSE, overwrite=TRUE, skip=TRUE)
        }
        SDMcurrent<-raster(paste(OutDirFinal, "/", crit,"/", sp, "/", sp, "_current.asc", sep=""), crs="+init=epsg:4326")
        SDMcurrent<-crop(SDMcurrent, myext)

        #calculate median, Q90 and Q10 for the species, save & zip
        Sys.sleep(0.1)
        print(paste(Sys.time(), "Median", sep=" "))
        Med<-calc(ascstack, fun=function(x) {quantile(x,probs = 0.5,na.rm=TRUE)})
        SuitChangeMed<-Med-SDMcurrent
        writeRaster(SuitChangeMed, filename=paste(SpOutDirFin,"/",sp,"_",year,"_MedianSuitChange.asc",sep=""), format="ascii", overwrite=TRUE)
        writeRaster(Med, filename=paste(SpOutDirFin,"/",sp,"_",year,"_MedianCrop.asc",sep=""), format="ascii", overwrite=TRUE)
        values(Med)<-ifelse(values(Med)<thresh, 0, values(Med))  #make all SDM raster values 0 if below threshold, 1 otherwise
        writeRaster(Med, filename=paste(SpOutDirFin,"/",sp,"_",year,"_MedianCropThresh.asc",sep=""), format="ascii", overwrite=TRUE)
        
        Sys.sleep(0.1)
        print(paste(Sys.time(), "Q90", sep=" "))
        Q90<-calc(ascstack,fun=function(x) {quantile(x,probs = 0.9,na.rm=TRUE)})
        SuitChangeMax<-Q90-SDMcurrent
        writeRaster(SuitChangeMax, filename=paste(SpOutDirFin,"/",sp,"_",year,"_Q90SuitChange.asc",sep=""), format="ascii", overwrite=TRUE)
        writeRaster(Q90, filename=paste(SpOutDirFin,"/",sp,"_",year,"_Q90Crop.asc",sep=""), format="ascii", overwrite=TRUE)
        values(Q90)<-ifelse(values(Q90)<thresh, 0, values(Q90))  #make all SDM raster values 0 if below threshold, 1 otherwise
        writeRaster(Q90, filename=paste(SpOutDirFin,"/",sp,"_",year,"_Q90CropThresh.asc",sep=""), format="ascii", overwrite=TRUE)
        
        Sys.sleep(0.1)
        print(paste(Sys.time(), "Q10", sep=" "))
        Q10<-calc(ascstack,fun=function(x) {quantile(x,probs = 0.1,na.rm=TRUE)})
        SuitChangeMin<-Q10-SDMcurrent
        writeRaster(SuitChangeMin, filename=paste(SpOutDirFin,"/",sp,"_",year,"_Q10SuitChange.asc",sep=""), format="ascii", overwrite=TRUE)
        writeRaster(Q10, filename=paste(SpOutDirFin,"/",sp,"_",year,"_Q10Crop.asc",sep=""), format="ascii", overwrite=TRUE)
        values(Q10)<-ifelse(values(Q10)<thresh, 0, values(Q10))  #make all SDM raster values 0 if below threshold, 1 otherwise
        writeRaster(Q10, filename=paste(SpOutDirFin,"/",sp,"_",year,"_Q10CropThresh.asc",sep=""), format="ascii", overwrite=TRUE)

        SDMstack<-stack(SDMcurrent,Q10,
                        Med,SuitChangeMed,
                        Q90,SuitChangeMax)
        titles<-c("Current suitability","future minimum suitability",
                  "future median suitability", "median (most likely) change in snake suitability",
                  "future maximum suitability", "highest potential increase in suitability")
        
        #remove unzipped ascis   
        #for (file in outfiles){
        #  file.remove(file)
        #}
        
        
        #print out summary of whats going on (keep track of progress: start)
        Sys.sleep(0.1)
        print(paste(Sys.time(),"done summarizing raw future projections for", sp,sep=" "))

        
        
        #######################################
        ############### Plotting ##############
        #######################################
        
        #print out summary of whats going on (keep track of progress: start)
        Sys.sleep(0.1)
        print(paste(Sys.time(),"plotting future sutability for", sp, "for",year,sep=" "))
        
        mybgxs<-crop(mybg, y=myext)
        
        png(filename=paste(PlotDir,"/",sp,"_",year,"_SDMFut.png",sep=""), units="in", width=10, height=12, res=500) #png smaller file size hat pdf
        par(mfrow = c(3,2), mar = c(1,1,1,1), oma=c(1,1,5,1), mgp=c(1.2,0.5,0))
        #mar=c(bottom, left, top, right); oma= outer margins; mpg = axis label locations c(3, 1, 0)
        
        for (p in c(1:6)){
          
          print(paste(Sys.time(),"plot",p ,sep=" "))
          
          plot(mybgxs, col=mycols1, ylab="latitude", xlab="longitude", box=FALSE, bty="l", cex.axis=0.8, legend=FALSE)
          
          if(p %in% c(1,2,3,5)){
            plot(SDMstack[[p]], col=my_col2, zlim= c(0,1), box=FALSE, bty="l", legend=TRUE,  add=TRUE)
          }else{ #plot change raters
            if(p %in% c(4,6)){
              plot(SDMstack[[p]], col=my_col5, zlim=c(-1,1), box=FALSE, bty="l", legend=TRUE,  add=TRUE)
            }}
        
          plot(worldmap[5], border="grey20", lwd= 0.5, add=TRUE,col=NA,reset=FALSE)
          
          for (mapunit in DatSum[is,"gen_sp_subtax"]){
            k<-which(DatSum$gen_sp_subtax==mapunit)
            q<-which(DatSum[is,"gen_sp_subtax"]==mapunit)
            yesCountry<-str_split(DatSum[k,"country"], "/")                   #split countries of occurrence into separate elements
            mycountries <- subset(worldmap, worldmap$GU_A3 %in% yesCountry[[1]])     #get subset of countries that species occurs in
            plot(mycountries[5], add=TRUE, col=NA, border="black",lwd= 1)
            
            species_shapeA <- subset(WHOpolies, WHOpolies$Scientific == gsub("_"," ",paste(mapunit)))
            
            if(length(species_shapeA)>=1){
              plot(species_shapeA[2], add=TRUE, border=mycols2[q], col=WHOcol, lwd= 1)
            }
            
            if(file.exists(paste(OCCDir,"/",mapunit,".csv",sep=""))){
              urecs<-read.table(paste(OCCDir,"/",mapunit,"_allRecs.csv",sep=""), sep=",", header=TRUE)
              #points(urecs$lat~urecs$long, col=mycols2[q], pch=1, cex=1)
            }
            if(file.exists(paste(SWDDir,"/",mapunit,".csv",sep=""))){
              urecs2<-read.table(paste(SWDDir,"/",mapunit,".csv",sep=""), sep=",", header=TRUE)
              #points(urecs2$lat~urecs2$long, col=mycols2[q], pch=3, cex=1)
            }
            
            if(file.exists(paste(OCCDir,"/",mapunit,".csv",sep=""))){
              if(length(urecs[,1])>=1){
                myx<-mean(urecs$long)
                myy<-min(urecs$lat)-2
                graphics::text(x=myx, y=myy, labels = gsub("_"," ",mapunit), cex = 1, col = mycols2[q])
              }
            }
            
          }
          
          legend(legend=paste(titles[p]),x="top", cex=1, bg="white",inset=c(0.02,0.02))
          
        }
        
        
        ############
        #plot title#
        ############
        
        mtext(paste("'",gsub("_"," ",spX),"'"," modelling unit",sep=""), side = 3, line = 1, outer = TRUE, cex= 1.5)
        
        dev.off()
        
        #print out summary of whats going on (keep track of progress: finish)
        Sys.sleep(0.1)
        print(paste(Sys.time(),"done plotting", sp,sep=" "))
       
      } 
      } 
    }
    ###################
  }
    ###################
} #end for species loop


############################################
############################################
# next: 8-VAP_CDCut ########################
############################################


